%% This file replicates the simulation results of the Chen and Szroeter (2012) test  (Table 3 and 4). For the results of the Bennett test see the R file  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Continous Outcomes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  a .2  b 0 n 250 
n=250;  
a=.2;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a .6  b 0 n 250
n=250;  
a=.6;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))


%%  a .2  b 1 n 250
n=250;  
a=.2;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a .6  b 1 n 250
n=250;  
a=.6;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))

%%  a .2  b 0 n 1000
n=1000;  
a=.2;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))



%%  a .6  b 0 n 1000
n=1000;  
a=.6;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))


%%  a .2  b 1 n 1000
n=1000;  
a=.2;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a .6  b 1 n 1000
n=1000;  
a=.6;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;
A=[-Inf;-1;0;1;Inf];

des=0;
Rc=zeros(M,1);
Rcp=zeros(M,1);
Rcp2=zeros(M,1);


P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=b1*D+b2*Z+U;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);
Rcp2(i)=chszpsim(Y,Z,D,P,K,B,A,cl);


end

res=[mean(Rc) mean(Rcp) mean(Rcp2)];
fprintf(' %1.4f  & %1.4f  & %1.4f  \\\\   \n', res(1,:))
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Binary Outcomes%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%  a .2  b 0 n 250
n=250;  
a=.2;
b2=0;
M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a  .6  b 0 n 250
n=250;  
a=.6;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))

%%  a  .2  b 1 n 250 
n=250;  
a=.2;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a  .6  b 1 n 250
n=250;  
a=.6;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))

%%  a  .2  b 0 n 1000
n=1000;  
a=.2;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))



%%  a  .6  b 0 n 1000
n=1000;  
a=.6;
b2=0;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))


%%  a  .2  b 1 n 1000
n=1000;  
a=.2;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))




%%  a  .6  b 1 n 1000
n=1000;  
a=.6;
b2=1;

M=1000;
b1=1;
sigma=1;
B=499;


cl=0.05;


des=1;
Rc=zeros(M,1);
Rcp=zeros(M,1);



P=2;
K=sqrt(n/(2*log(log(n))));




parfor i=1:M
rng(i);
Z=rand(n,1)>0.5;
UV=mvnrnd([0;0],[sigma 0.5;0.5 sigma],n);
U=UV(:,1);
V=UV(:,2);
D=(a*Z+V)>0;
Y=(b1*D+b2*Z+U)>0;
Y=double(Y);
D=double(D);
Z=double(Z);

A1=min(Y):(max(Y)-min(Y))/2:max(Y)+0.00001;
Rc(i)=chszsim(Y,Z,D,P,K,des,cl);
Rcp(i)=chszpsim(Y,Z,D,P,K,B,A1,cl);



end

 res=[mean(Rc) mean(Rcp)];
fprintf(' %1.4f  & %1.4f  \\\\   \n', res(1,:))